/*
 * OREKIT-X
 * Copyright 2002-2008 CS Communication & Systemes
 * 
 * Parts of this software package have been licensed to CS
 * Communication & Systemes (CS) under one or more contributor license
 * agreements.  See the NOTICE file distributed with this work for
 * additional information.
 *  
 * This is an experimental copy of OREKIT from www.orekit.org.
 * Please use the original OREKIT from orekit.org for normal work
 * unrelated to this research project.
 * 
 * Licensed under the Apache License, Version 2.0 (the "License"); you
 * may not use this file except in compliance with the License.  You
 * may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
 * implied.  See the License for the specific language governing
 * permissions and limitations under the License.
 */
package ore;

import ore.errors.OrekitException;
import ore.time.AbsoluteDate;
import ore.time.TimeStamped;

import org.apache.commons.math.geometry.Vector3D;

/**
 * This class handles orbital parameters without date.
 * 
 * The aim of this class is to separate the orbital parameters from the date
 * for cases where dates are managed elsewhere. This occurs for example during
 * numerical integration and interpolation because date is the free parameter
 * whereas the orbital parameters are bound to either differential or
 * interpolation equations.
 * 
 * For user convenience, both the cartesian and the equinoctial elements
 * are provided by this class, regardless of the canonical representation
 * implemented in the derived class (which may be classical keplerian
 * elements for example).
 * 
 * The parameters are defined in a frame specified by the user. It is important
 * to make sure this frame is consistent: it probably is inertial and centered
 * on the central body. This information is used for example by some
 * force models.
 * 
 * An instance of this class is immutable.
 * 
 * @author Luc Maisonobe
 * @author Guylaine Prat
 * @author Fabien Maussion
 * @author V&eacute;ronique Pommier-Maurussane
 */
public abstract class Orbit 
    extends Object
    implements TimeStamped, 
               java.io.Serializable 
{

    public final Frame frame;

    public final AbsoluteDate date;

    public final double mu;

    private volatile PositionVelocity pvCoordinates;


    /** 
     * @param frame Frame in which the parameters are defined
     * @param date Date for which the orbital parameters are defined
     * @param mu Central attraction coefficient (m^3/s^2)
     */
    protected Orbit(Frame frame, AbsoluteDate date, double mu) {
        super();
        this.frame =  frame;
        this.date = date;
        this.mu = mu;
    }
    /**
     * Define an orbit from cartesian parameters.
     * @param pvCoordinates the position and velocity in the inertial frame
     * @param frame the frame in which the {@link PositionVelocity} are defined
     * @param date date of the orbital parameters
     * @param mu central attraction coefficient (m^3/s^2)
     */
    protected Orbit(PositionVelocity pvCoordinates, Frame frame,
                    AbsoluteDate date, double mu)
    {
        super();
        this.frame = frame;
        this.date = date;
        this.mu = mu;
        this.pvCoordinates = pvCoordinates;
    }
    protected Orbit(Orbit orbit){
        super();
        this.frame = orbit.frame;
        this.date = orbit.date;
        this.mu = orbit.mu;
    }


    public final Frame getFrame() {
        return this.frame;
    }
    public final double getMu() {
        return this.mu;
    }
    public final AbsoluteDate getDate() {
        return this.date;
    }
    /** 
     * @return Semi-major axis (m)
     */
    public abstract double getA();
    /** 
     * @return First component of the equinoctial eccentricity vector
     */
    public abstract double getEquinoctialEx();
    /** 
     * @return Second component of the equinoctial eccentricity vector
     */
    public abstract double getEquinoctialEy();
    /** 
     * @return First component of the inclination vector
     */
    public abstract double getHx();
    /** 
     * @return Second component of the inclination vector
     */
    public abstract double getHy();
    /** 
     * @return Eccentric latitude argument (rad)
     */
    public abstract double getLE();
    /** 
     * @return True latitude argument (rad)
     */
    public abstract double getLv();
    /** 
     * @return Mean latitude argument (rad)
     */
    public abstract double getLM();
    /** 
     * @return Eccentricity
     */
    public abstract double getE();
    /** 
     * @return Inclination (rad)
     */
    public abstract double getI();
    /** 
     * The keplerian period is computed directly from semi major axis
     * and central acceleration constant.
     * @return Keplerian period in seconds
     */
    public double getKeplerianPeriod() {
        double a = this.getA();
        return (2.0 * Math.PI * a * Math.sqrt(a / this.mu));
    }
    /** 
     * The keplerian mean motion is computed directly from semi major
     * axis and central acceleration constant.
     * @return Keplerian mean motion in radians per second
     */
    public double getKeplerianMeanMotion() {
        double a = this.getA();
        return Math.sqrt(this.mu / a) / a;
    }
    /** 
     * @param outputFrame frame in which the position/velocity coordinates shall be computed
     * @return pvCoordinates in the specified output frame
     * @exception OrekitException if transformation between frames cannot be computed
     * @see #getPositionVelocity()
     */
    public PositionVelocity getPositionVelocity(final Frame outputFrame)
        throws OrekitException
    {
        if (this.pvCoordinates == null) {
            this.initPositionVelocity();
        }

        /* 
         * If output frame requested is the same as definition frame,
         * PV coordinates are returned directly
         */
        if (outputFrame == frame)
            return pvCoordinates;
        else {
            /* 
             * PV coordinates are transformed to output frame
             */
            Transform t = frame.getTransformTo(outputFrame, date);
            return t.transformPositionVelocity(pvCoordinates);
        }
    }
    /** 
     * @return pvCoordinates in the orbital definition frame
     * @see #getPositionVelocity(Frame)
     */
    public PositionVelocity getPositionVelocity() {
        if (this.pvCoordinates == null) {
            this.initPositionVelocity();
        }
        return this.pvCoordinates;
    }
    /** 
     * Initialize the position/velocity coordinates.
     */
    private void initPositionVelocity() {
        /*
         * Equinoctial parameters
         */
        double a  = getA();
        double ex = getEquinoctialEx();
        double ey = getEquinoctialEy();
        double hx = getHx();
        double hy = getHy();
        double lE = getLE();
        /*
         * Inclination related intermediate parameters
         */
        double hx2   = hx * hx;
        double hy2   = hy * hy;
        double factH = 1. / (1 + hx2 + hy2);
        /* 
         * Reference axes defining the orbital plane
         */
        double ux = (1 + hx2 - hy2) * factH;
        double uy =  2 * hx * hy * factH;
        double uz = -2 * hy * factH;

        double vx = uy;
        double vy = (1 - hx2 + hy2) * factH;
        double vz =  2 * hx * factH;
        /*
         * Eccentricity-related intermediate parameters
         */
        double exey = ex * ey;
        double ex2  = ex * ex;
        double ey2  = ey * ey;
        double e2   = ex2 + ey2;
        double eta  = 1 + Math.sqrt(1 - e2);
        double beta = 1. / eta;
        /*
         * Eccentric latitude argument
         */
        double cLe    = Math.cos(lE);
        double sLe    = Math.sin(lE);
        double exCeyS = ex * cLe + ey * sLe;
        /*
         * Coordinates of position and velocity in the orbital plane
         */
        double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - ex);
        double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - ey);

        double factor = Math.sqrt(this.mu / a) / (1 - exCeyS);
        double xdot   = factor * (-sLe + beta * ey * exCeyS);
        double ydot   = factor * ( cLe - beta * ex * exCeyS);
        /*
         */
        Vector3D position =
            new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);

        Vector3D velocity = 
            new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);

        this.pvCoordinates = new PositionVelocity(position, velocity);
    }
    public final int compareTo(TimeStamped stamped){
        return this.date.compareTo(stamped);
    }
}
